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AN ADAPTIVE FINITE ELEMENT METHOD FOR HIGH STEED COMPRESSIBLE FLOW 


i ' 


R. Lohner. K. Morgan ar.d O.C. Zienkiewicz 

Institute for Numerical Methods art Engineering, University College of Swansea, 
Singleton Park, Swansea SA2 8FP, U.K. 

Our aia is to solve large 2/3-D compressible fluid flow problems employing 
finite elements. Other numerical techniques, like finite differerce or finite 
volume methods have reached a high degree of sophistication [1,2], but it is ex- 
pected that the finite element method will make significant contributions due to its 
geometrical flexibility, a factor that is of importance for industrial applications. 

The development has so far gone through the following stages: 

(a) Definition of the basic algorithm: in order to be competitive, end at the same 
time to model correctly the physical properries of hyperboi’c equations, it was 
decided to employ explicit time-marching schemes. As the straightforward Galerkin 
method is suboptimal (3], a modified or ’upwinded' Taylor Galerkin-type (4) pro- 
cedure was implemented. The complete description may be found in (5], together 
with numerical examples. This one-step algorithm has now been superseded by an 
equivalent two-step method, the description of which may be found in [ 6 ] . This 
two-step scheme is 2-3 times faster than the one-step scheme on scalar machines 
(for two-dimensional problems) and has been vectorized in order to realise the full 
power of modern supercomputers. 

(b) Domain splitting: it is well known that solutions of high speed compressible 

flow problems exhibit narrow regions of rapid change (e.g. shocks) which are em- 
bedded in lerger regior.. where the solution is smooth. Accordingly, large vari- 
ations in element s’zo are expected in typical discretizations. However, the 
small elements might then require that a correspondingly small global timestep 

co Id be employed in larger elements. The remedy adopted here, and described in 
det.il in ( 7 j , is to split the domain into regions in which different timestep- 
sizes can be used. The domain subdivision is performed completely automatically 
by the computer code at prescribed time. intervals, and allows a time-accurate dev- 
elopment of the unsteady solution. 

(c) Adaptive mesh refinement : in general, an analyst will have no a priori know- 

ledge of the location of those areas of the domain where more (i.e. smaller) el- 
ements should be employed. Therefore, usually, much more elements than necessary 
will be employed, leading to an inefficient overall procedure. An ideal comput- 
ational algorithm would require the ability to refine the mesh where necessary as 
the solution proceeds. The geometric flexibility of the linear triangular el- 
ement makes it ideally suited for refinement processes of this type. We adopted a 
posteriori methods [8), as they seem at present more economical, and for th» ame 
reason also did not implement hierarchical techniques |9], but the more classic 
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enrichment of addi»( sore elwnts. At present we are only considering steady 
state problems. the generalization to transient problems being an obvious eatensior. 
This wans, that Che timestepping scheme is utilized as a relaxation procedure. 

After a given number of timesteps the solution doaain is analysed, and aore el— 
eaents are added where necessary. Cenerally speaking, there exist three poss- 
ibilities for approxinately determining the error 


where u denotes the exac and u the discrete solution: 


i ) Comparison with h order schemes: The significant derivatives of the 

partial differential equation (PDE) under consideration are evaluated twice, using 
in each case difference schemes of different order [10, II, 12). By determining 
the discrepancy of both approximations an estimation of the error can then be ob- 
tained. The problem with this kind of approach is that it fails near boundaries 
and at singularities or boundary layers (which are common in fluid dynamic prob- 
lems). At the same time it is not extendable to FEMs, which operate on an element 


ii) Determination of the relative importance of further degrees of freedom: 


Further degrees of freedom are introduced on an element by element oasis, and the 
relative importance of these further degrees gives an error estimate (13, 14, IS]. 

The problem with this kind of approach is that it is relatively expensive in CPU- 
time requirements, so that for transient nroblems a considerable percentage of run- 
time will be spent on error estimation. 

iii) Use of error norms: Here the classic theoretical error estimates are employed 

locally [ 16.17], Thus, no further degrees of freedom are introduced and only 
first or second derivatives need to be evaluated. Our experience indicates chat 
this type of error indicator works satisfactorily, and, as it is very economical, 
it is regarded as a good algorithm for transient problems as well. For elliptic 
problems the appropriate error norms appear naturally whereas for hyperbolic prob- 
lems the theory is far from complete. Nevertheless one can assume 


1 “ * < c h* k |u| r 


where h is a representative element length. Using the L.-norro (k-0) yields 


1“ ~ < c h £ |u| f . 


The aim of any refinement is to obiain a reduction of errors according to some 
criterion, e.g. at a certain point, surface or evenly throughout the field. Par- 
ticularly for hyperbolic problems the error at one point may influence the accuracy 
ol the solution in the whole field (e.g. the root of an expansion fan), so that an 
even distribution of errors seems to be the only possible practical choice. 
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Therefore, at each refinement level all elements satisfying 


h^ju| > a max h*' |u|. W) 

e ■ 1 .nelem 

are refined. Sinre the exact solution u is unknown, the practical requirement 
becomes 

h^|u**| > a cun h^ | «a** I m ■ 

1 e*l, nelem 

Only the cases 1=1 or 1-2 appear to be of practical interest, and both have been 
studied (see examples). For the case 1-2 the first derivatives of a are evaluated 
inside the elements, and hereafter the nodal values for the second derivatives are 
recovered variationally as follows: 

| N l dV - - | N 1 M V dV , (6) 

where M* 1 is constant and u,* 1 is defined on an element basis. It has been found 

x 

that a-values of the order 


a - 0.6 - 0.9 (7) 

yield the most effective refinement strategy. This is in contrast to (13), where 
tne factor a - 0. 1 was reported as optimal. A possible explanation for the dis- 
crepancy of these values may be found in the nature of the PDFs treated in both 
cases: whereas here the PDEs are hyperbolic - and this means th-t small disturb- 
ances propagate far into the field - , in [12) the effective solution of elliptic 
PDEs was pursued - and this means that small disturbances decay rapidly. 


Resul ts 

(a) Superson ic flow past a wedge: the successive stages of the domain discretiz- 
ation as well as the solution obtained are shown in figire 1 . In this case the 

mesh was enriched according to equation (5) with 1-1 end a-0.6. 

(b) Prandtl-Meyer expansion fan: the problem statement, as well as the successive 

stages of the domain discretization and the corresponding solutions are depicted 
in figure 2. The improvement in solution quality is readily seen. In this 

case the mesh was enriched according to equation (5) wi.\ £-2 and a-0.8. 
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